Estimation of Landsat-like daily evapotranspiration for crop water consumption monitoring using TSEB model and data fusion

Evapotranspiration (ET) plays an essential role in agricultural water resource management. Understanding regional agricultural water consumption characteristics can be improved by predicting ET using remote sensing. However, due to the lack of high-resolution images on clear-sky days or the limitation of ET reconstruction on cloudy-sky days, it remains challenging to continuously derive ET at the field scale. In this study, the Landsat and MODIS data were initially fused to obtain the Landsat-like vegetation index and land surface temperature on clear-sky days. Then the two-source energy balance (TSEB) model was applied to calculate the daily ET during the clear-sky. A canopy resistance-based gap-filling method was involved in reconstructing regional ET on cloudy days while considering different environmental factors. The estimations were validated by automatic weather system data (AWS) and eddy covariance (EC) measurements in Guantao County. The results demonstrated that the proposed scheme performed well in estimating cropland ET, with an RMSE of 0.86 mm·d−1 and an R2 of 0.65, and the NSE and PBias were 0.61 and -0.29%, respectively. The crop water consumption analysis revealed that the daily ET of winter wheat peaked during the maturation stage. Nevertheless, summer maize water consumption peaked in the middle of the growing season in this area. The temperature during the early development stage and the soil moisture in the mid and late growth stages had the greatest impact on the ET of winter wheat. During the entire growing period, soil moisture had the largest effect on the ET of summer maize. The findings showed that the TSEB model can be effectively applied to field-scale water consumption monitoring in North China through MODIS and Landsat data fusion and ET temporal reconstruction considering environmental factors.


Introduction
Evapotranspiration (ET) is an essential medium of the interaction between the surface and the atmosphere [1,2]. It plays a crucial role in agricultural systems' energy and water balance. This study intended to generate daily Landsat-like ET using the TSEB model for crop water consumption through MODIS and Landsat data fusion in North China. Section 2 introduces the research area, remote senisng data and station measurements. Then the data fusion scheme and ET temporal reconstruction method for TSEB were proposed. In Section 3, the performance of ET estimates was systematically evaluated with surface energy observations and daily ET measurements from AWS and EC stations. The crop water consumption was analyzed according to the Landsat-like ET. Section 4 discussed the uncertainty of the TSEB model and analyzed the impact of environmental factors on crop water consumption in different growing seasons. In Section 5, we conclude the achievements for crop water consumption monitoring in Guantao County, North China.

Study area
The North China Plain (N35.02˚-N40.41˚, E113.18˚-E119.82˚) is one of the three major plains in China, with the largest population and acts as main grain-producing region in China. Guantao County is located in the North China Plain, with a semiarid and semihumid mainland monsoon climate. The yearly mean temperature, annual potential evaporation, and precipitation are approximately 14˚C, 1100.65 mm, and 548.7 mm, respectively.
Due to the large amount of agricultural water consumption and insufficient precipitation, crop irrigation water is facing a significant shortage. Therefore, a large amount of groundwater is exploited in this area to supplement irrigation. The main planting form in this area is the rotation of summer maize and winter wheat. The winter wheat mainly grows from October to May, and the summer maize grows from June to September. The land-use types for the study area are shown in Fig 1, and cropland is the dominant land type in this area. The red star in Fig 1 is the automatic weather station (AWS) and eddy covariance (EC) station.

EC and AWS measurements
The EC and AWS stations' longitudes and latitudes are 115.1 E, 36.5 N, respectively. The ground-based measurements from 2008 to 2010 were used to validate the results of fused Landsat-like LST and simulated daily ET [26,27]. EdiRe software was adopted to deal with the raw EC observations, including spike detection, lag correction of H 2 O/CO 2 relative to the vertical wind component, sonic virtual temperature correction, coordinating rotation, correction for density fluctuation, and frequency response correction. The final released EC measurements provided half-hourly sensible heat flux (H) and latent heat flux (LE).
The AWS measured the surface net radiation (Rn), wind speed, air temperature, and humidity. The soil heat flow plate was buried underground to measure the soil heat flux (G), mean soil temperature, and soil moisture. The surface temperature from AWS was measured by a thermal infrared thermometer. The AWS observations were averaged to half-hour frequency to be consistent with the EC data.
The Bowen ratio correction method was used to solve the energy non-closure, which assumes an invariable specific value between latent heat and sensible heat in the whole closure correction process [28]. Last, continuous daily ET measurements can be obtained by applying the mean diurnal variation (MDV) method to the EC latent heat gap-filling [29].

The Satellite and GLDAS data
The LST data from the MOD11A1 product at 1 km spatial resolution and land surface reflectance data from the MOD09GA at 500 m spatial resolution were applied to extract fusion parameters, i.e., the surface radiometric temperature and normalized difference vegetation index (NDVI).
The spatial resolution of Landsat 5 is 30 m in optical bands and 120 m in the thermal infrared band. The data were atmospherically corrected by the Pixel Information Expert (PIE) remote sensing image processing platform to derive the surface reflectance. The surface albedo α was computed according to Liang et al. [30]. The Landsat red and near-infrared bands' reflectance was used to calculate NDVI. In this research, the surface radiometric temperature T rad of Landsat 5 was computed by band six based on the single-channel algorithm put forward by Sobrino et al. [31]. The atmospheric correction parameters can be achieved on the NASA official website from the Atmospheric Correction Parameter Calculator.
The Global Land Data Assimilation System (GLDAS) contains various global-scale meteorological elements. The data used in this article included surface air pressure, specific humidity, soil moisture content (0-10 cm), downward shortwave radiation flux, air temperature, and wind velocity. We resampled all the data to a 30 m spatial resolution by bilinear interpolation. Although GLDAS has a coarse resolution, it can still provide meteorological element values of several pixels in the study area, which is more spatially representative than the information of one meteorological station.

Methods
This research applied the spatial and temporal adaptive reflectance fusion model (STARFM) to obtain the NDVI and LST from clear days' Landsat and MODIS data [32]. The acquisition of fused NDVI and LST using the same pixel searching processes. Then, the clear days' ET and canopy resistances were calculated and inversed according to the TSEB and Penman-Monteith (PM) models, respectively [33]. The canopy resistances of the cloudy days were reconstructed to obtain continuous crop ET at the field scale by the daily environmental factors and the clear days' canopy resistances. Fig 2 shows the overall method flow chart of this research.
2.4.1 Data fusion model. The relatively high spatial resolution remote sensing image, e.g., Landsat 5, has a long revisit cycle. The cloud contamination of images resulted in no observation records available within the study area during the growing season. The MODIS satellite sensor can obtain overpass data every day, but the spatial resolution is low. Gao et al. [32] proposed STARFM to fuse surface reflectance from Landsat and MODIS. The basic idea of the model is to combine the advantages of Landsat's high spatial resolution and MODIS's high temporal resolution. This research fused the MODIS and Landsat NDVI and LST from 2008 to 2010 to obtain clear day data at a Landsat-like spatial resolution.
2.4.2 TSEB model. The two-source energy balance (TSEB) model has been applied worldwide [34]. It uses the net radiation Rn, sensible heat flux H, and latent heat flux LE as energy

PLOS ONE
Estimation of Landsat-like daily evapotranspiration for crop water consumption monitoring using TSEB PLOS ONE | https://doi.org/10.1371/journal.pone.0267811 May 19, 2022 components. The relationship between them is as follows: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2 cosðy z q ÞÞ� ð3Þ where R nc and R ns are the net radiation over the canopy and soil (W�m −2 ); H c and H s are the sensible heat flux over the canopy and soil, respectively; LE c and LE s are the canopy and soil radiation, respectively; and G is the soil heat flux. It is supposed to be an invariable ratio of R ns [12], and subscripts c and s are defined as vegetation and soil, respectively. S d is the downwelling shortwave radiation (W�m −2 ); ε is the land surface emissivity; ε a is the atmospheric emissivity computed by the air temperature T a (K) and water vapor pressure [35]. In Eq (3), LAI is the leaf area index and has a linear relationship with NDVI [12]. k is an extinction coefficient for LAI related to the canopy structure [36]. θ z is the solar zenith angle. c g is assumed to be a constant ratio equal to 0.35. H c and H s are computed by Eqs (7) and (8), ρ is the air density (kg�m −3 ), C p is the heat capacity of air under constant pressure(J�kg −1 �K −1 ), and T c and T s are the radiometric temperatures from the canopy and the soil surface, respectively. r a,c is the aerodynamic resistance to heat transfer between the canopy and the reference height, and r a,s is the aerodynamic resistance over the soil surface according to Morillas et al. [37]: ð9Þ f c is estimated by the fractional vegetation cover considering the sensor view angle. T c and T s are the canopy and soil radiometric temperatures, respectively. The TSEB uses the Priestley-Taylor (P-T) method (Equations (4) and (7)) to calculate the initial canopy temperature T ci : where f g is the fraction of the green vegetation, α pt has an initial value of 1.26. γ is the psychrometric constant of~0.066 (kpa�K −1 ), and Δ is the slope of the saturation vapor pressure versus the air temperature curve. H c and H s can be calculated by combining Eqs (7) and (8) and (9) with the initial T c , and LE s are computed by Eqs (1)- (5). If the value of LE s is negative, an iterative process is adopted to produce a new estimate of T c by the decline of α pt .
The evaporative fraction (EF) method was used to extend the instantaneous evapotranspiration to the daily scale. It assumes that the evaporative fraction, i.e., LE/(Rn-G) is a constant value during the daytime. L is the latent heat of vaporization (MJ�m -2 �mm -1 ). The daily net radiation (R n24 , MJ�m -2 �d -1 ) was computed by the solar radiation from the GLDAS data.
G 24 is defined as the daily average soil heat flux and is approximately equal to zero. LE 24 is the daily latent heat flux.  [38]. m(Tmin) and m(VPD) were used to limit potential stomatal conductance by minimum air temperatures (Tmin) and reduce the potential stomatal conductance when VPD is high enough to inhibit photosynthesis, respectively. Therefore, by introducing the canopy resistances of cloudy days into the TSEB model, the ET on cloudy days can be obtained. Table 1 summarizes the parameters used in this study.

Landsat-like fused-LST
The LST fused by the STARFM model was compared with the ground observations. The determination coefficient (R 2 ) between fused-LST and the surface temperatures measured by AWS shown in Fig 3 was 0.95. The percentage bias (PBias) and the root mean square error (RMSE) were −0.33% and 2.55 K, respectively. The result indicated that although fused-LST cannot entirely reflect the actual temperature of the surface, the fused-LST obtained by STARFM has reasonable accuracy.

Daily net radiation
The daily net radiation flux has a considerable influence on the estimation of daily ET. The TSEB model used daily net radiation (Rn) as a primary boundary condition. The daily net radiation computed from the GLDAS solar radiation and surface albedo was compared with the AWS measurements. Fig 4 shows [39], the accuracy was regarded as reasonable if the NSE and R 2 were more significant than 0.5 and 0.6, respectively.

Validation of daily ET
The daily ET computed from TSEB was compared with the EC measurements from 2008 to 2010. As shown in Fig 5(a) and Table 2, the TSEB shows a good performance of daily ET compared to the land truth from EC observations. As a whole, the R 2 , NSE, and RMSE of the estimated daily ET at the field scale were 0.65, 0.61, and 0.86 mm, respectively, and had no The change characteristics of monthly ET obtained by TSEB and EC observations are demonstrated in Fig 5(b). The values of the NSE, R 2 , and PBias on the monthly scale became larger than those on the daily scale. These validation results showed that using the proposed method in this research can effectively estimate continuous ET at the field scale.

Crop water consumption
As demonstrated in Fig 6, the water consumption of summer maize and winter wheat during the whole growth period showed noticeable monthly changes. The middle growth period of summer maize was from July to August and was the most vigorous stage. The transpiration of leaves reached a maximum during this period. However, the maturation growth stage of winter wheat in June had the highest ET. Fig 6 shows that in the early period of growth, the average water consumption of summer maize was greater than that of winter wheat. Nevertheless, the average water consumption of summer maize in the middle and late growth stages was less than that of winter wheat. Table 3 shows the water consumption of summer maize and winter wheat at different growth stages.  Table 3 shows that the total water consumption of winter wheat was less than that of summer maize throughout the growth stages. Winter wheat had similar water consumption in different growth stages from 2008 to 2009 and 2009 to 2010. The total water consumption for the two periods was 267.85 mm and 291.36 mm, respectively. However, the water consumption of summer maize varied significantly in different growth stages in 2009 and 2010. The water consumption of summer maize in 2009 was substantially greater than that in 2010; with a total water consumption for the two periods of 367.04 mm and 300.67 mm, respectively. This increase in consumption is generally due to greater evapotranspiration caused by more irrigation of cropland. was much higher than that in August 2010. The reason might be that irrigation or rainfall increases the soil moisture, and soil evaporation and vegetation transpiration become more significant. The difference in ET between irrigated and nonirrigated croplands will increase. Overall, the water consumption of summer maize is greater than that of winter wheat. In the North China Plain, rainfall is mainly concentrated during summer, and the winter wheat is more dependent on irrigation. There were also differences in water consumption between different years for both winter wheat and summer maize, which may be related to the crop planting density or irrigation method.

Discussion
High-resolution ET data can provide more details for the water consumption of plants. The application of the TSEB model in field-scale crop water consumption monitoring in North China still needs careful investigation. This study used data fusion and the TSEB model, which considers environmental factors, to estimate daily field-scale ET in the North China Plain.
The LST plays a significant role in ET estimation, but the data fusion model has uncertainty. The validation result in Fig 3 shows that the fused-LST obtained in this paper has reasonable accuracy. This is the premise that the subsequent model can perform well. As shown in Fig 4, the daily net radiation is overestimated. The reason for the overestimation of daily net radiation is generally the uncertainty of meteorological data. According to the TSEB model, the daily net radiation is inversely proportional to albedo. Although there may be a gap between meteorological data and actual values, from the results, the accuracy of daily net radiation is still within an acceptable range. The overestimation of daily net radiation also leads to the overestimation of overall ET in Table 2. The overall validation shows that the accuracy of the monthly statistical results is significantly higher than that of the daily statistical results, which explains why the estimation error of the scheme proposed by the article tends to compensate over time. This is mainly due to the random errors associated with daily evapotranspiration cancels each other out after time accumulation. Compared with previous ET research, which can only be carried out on days with clear-sky, the temporal reconstruction method proposed in this paper improved the scenario in which effective data cannot be obtained on cloudy-sky days. Through the data fusion model, the LST on cloudy-sky days was obtained, and then the ET on cloudy-sky days was estimated. When the effective LST cannot be obtained by using the data fusion model, the ET on cloudy-sky days was obtained by introducing the canopy resistance on cloudy-sky days into the TSEB model. This greatly improved the temporal resolution of ET research. The TSEB model that considers environmental factors can better analyze the variation law of the main environmental factors affecting evapotranspiration in the crop growth stages. The results indicated that the proposed scheme could conduct accurate ET estimation at a high spatiotemporal resolution at the field scale. However, reconstruction of canopy resistance on cloudy days using clear days calculation may affect the results because the growth of crops does not fit this general linear relationship. A time interval of at least five days is indispensable for accurate daily ET reconstruction, according to Alfieri et al. (2017). Obviously, errors may occur if the gap size is more extensive [40]. Irrigation also affects evapotranspiration, increasing the irrigation amount, soil water content, and evapotranspiration. Many previous studies have shown that ET is connected to crop species and may be influenced by environmental factors. These environmental factors are complex, and the common ones are solar radiation, air humidity, air temperature, and soil moisture [10,41].
The TSEB model includes a variety of exchange patterns between canopy and soil latent heat flux due to vegetation stress rather than directly analyzing the impact of environmental factors. The feedback of the TSEB model has a difference in environmental stress in semihumid and semiarid regions. Jarvis's type equations [42,43] were chosen to evaluate the model sensitivity and analyze the influence of environmental factors on crop growth. It contains four environmental factors, i.e., solar radiation (f sr ) soil moisture (f sm ), air humidity (f hu ), and air temperature (f ta ). The f sr factor describes the influence of solar radiation and implies the energy stress for heat flux [43]. It is calculated from downward shortwave radiation. The f sm factor considers the effect of soil moisture on canopy transpiration. The factor f hu is related to water stress conditions and addresses the impact of vapor pressure difference (VPD) in the atmosphere. The f ta factor represents the temperature constraint on canopy resistance. It describes the correlation between the air temperature at the reference height level over the canopy and the optimal air temperature (e.g., 298 K) for crop growth. Theoretically, the values of all environmental factors are between 0-1.
The mean environmental factor values for the summer maize and winter wheat throughout the entire growth stage are shown in Table 4. Winter wheat and summer maize's corresponding mean environmental factors such as f sm , f hu , and f ta differed significantly. This result indicated that the winter wheat was more susceptible to f sm , f hu , and f ta when the growing environment was drier and colder than that of summer maize. In contrast, the f sr for winter wheat and summer maize did not significantly differ. Soil moisture was the primary influencing factor for summer maize and winter wheat throughout the whole growing season. Tables 5 and 6 show the values of the environmental factors for winter wheat and summer maize during the three growing periods. The f sm and f hu had similar matters in the early growth stage of summer maize (sowing stage). Summer maize was mainly affected by soil moisture and air humidity, most likely due to irrigation; this suggested that both factors affect water vapor transfer from the canopy to the atmosphere. The air temperature was appropriate for the middle period of the summer maize growth (seeding and flowering stages). Solar radiation gradually replaced air humidity as the main factor affecting summer maize, and only soil moisture became a significant influencing factor during the maturation stage of growth. Table 6 shows that in the early period of winter wheat (sowing stage), the air temperature became the major influencing factor in this period. During the middle growth period for the winter wheat (wintering, reviving, and flowering stages), the impact of air temperature decreased because the weather warmed. During the maturation growth of winter wheat, soil moisture and air humidity were identical, which jointly became the main influencing factor for this period. The solar radiation factor was relatively high throughout the summer maize and winter wheat, indicating that solar radiation energy had a weak effect on crop evapotranspiration. Although solar radiation played a key role, the soil moisture data were not introduced to estimate the cropland evaporation. The spatial resolution of remote sensing-based soil moisture products does not meet the requirements of estimated evapotranspiration.

Conclusions
Accurate and continuous estimation of cropland ET is essential for agricultural irrigation scheduling, a necessary measure for the rational utilization of water resources, monitoring droughts, and aiding in scarce water supply management. It is critical in the arid and semiarid areas of North China, where water resources are limited.
This study generated Landsat-like daily ET using data fusion combined with the TSEB model for crop water consumption in Guantao County, North China. Through the fusion of MODIS and Landsat data, the high-resolution surface temperature and vegetation index data on clear days were obtained, and then the ET on cloudy days was obtained using the ET reconstruction method considering environmental factors. The validation demonstrated that this ET computation scheme performed well in crop water consumption monitoring. The overall R 2 is 0.65, NSE is 0.61, RMSE is 0.86 mm�d −1 , and PBias is -0.29%. The analysis of crop water consumption showed that the water consumption of winter wheat and summer maize peaked during the mature stage and the middle of the growing season, respectively. Among them, for winter wheat, the largest factor affecting ET in the early development stage was temperature, while in the middle and late growth stages, soil moisture had the greatest impact on ET. For summer maize, soil moisture played a decisive role in the change in ET throughout the growth period. This study demonstrates the feasibility of the TSEB combined with a data fusion model to generate daily ET for crop water consumption monitoring with a high spatial resolution (30 m). This will be a good and feasible practice for the rational allocation of water resources in the North China Plain and other similar agricultural areas globally.